function IdSet=OneSide1_ind(u01grid, u11grid, u21grid,u02grid, u12grid, u22grid,...
                            equalorder, ...
                            Ugridreduced, sgreduced,...
                            id_Ugridreduced,...
                            s_id_gridreduced, sU_id_gridreduced,...
                            n_U_sets, Tcoord, indices_pairs,...
                            numeqconstraints,coordrowseq, fillAeq, beq)
%% Loop over parameter values 
exists=zeros(sgreduced,1);
eps=10^(-6);
for g=1:sgreduced
% Set indices to keep track of equal elements     
idU1=id_Ugridreduced{1}(g,:);
idU2=id_Ugridreduced{2}(g,:);
idU3=id_Ugridreduced{3}(g,:);
idU1_t=id_Ugridreduced{1}(g,:).';
idU2_t=id_Ugridreduced{2}(g,:).';
idU3_t=id_Ugridreduced{3}(g,:).';
% Relevant sizes of the vector of unknowns
s1=s_id_gridreduced(g,1);
s2=s_id_gridreduced(g,2);
%s3=s_id_gridreduced(g,3);
sU=sU_id_gridreduced(g); 

%%%%%%%%%%%%%%%%%%%%%% Equality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Coordcolumnseq
id_01=[changecoord(s1,s2,idU1((1:4)),idU2(4),idU3(1))...  
       changecoord(s1,s2,idU1(4), idU2((1:3)),idU3(1))...
       changecoord(s1,s2,idU1((1:4)),idU2(4),idU3(3))... 
       changecoord(s1,s2,idU1(4), idU2((1:3)),idU3(3))...
       changecoord(s1,s2,idU1((1:4)),idU2(4),idU3(4))... 
       changecoord(s1,s2,idU1(4), idU2((1:3)),idU3(4))...
       changecoord(s1,s2,idU1(1:4),idU2(1), idU3(2)) ...
       changecoord(s1,s2,idU1(1:4),idU2(2), idU3(2)) ...
       changecoord(s1,s2,idU1(1:4),idU2(3), idU3(2)) ...
       changecoord(s1,s2,idU1(1:4),idU2(4), idU3(2)) ...
       changecoord(s1,s2,idU1(3),idU2(3),idU3(1))].'; %ones (columns)


coordcolumnseq= [changecoord(s1,s2,idU1(1),idU2(3),idU3(3)) changecoord(s1,s2,idU1(3),idU2(3),idU3(3)) changecoord(s1,s2,idU1(1),idU2(3),idU3(1)) ... %observational equivalence
                 changecoord(s1,s2,idU1(3),idU2(3),idU3(3)) changecoord(s1,s2,idU1(3),idU2(1),idU3(3))...
                 changecoord(s1,s2,idU1(1),idU2(1),idU3(1)) ... 
                 changecoord(s1,s2,idU1(2),idU2(3),idU3(4)) changecoord(s1,s2,idU1(3),idU2(3),idU3(4)) changecoord(s1,s2,idU1(2),idU2(3),idU3(1)) ... %observational equivalence
                 changecoord(s1,s2,idU1(3),idU2(3),idU3(4)) changecoord(s1,s2,idU1(3),idU2(2),idU3(4))...
                 changecoord(s1,s2,idU1(2),idU2(2),idU3(1)) ...  
                 id_01.']; %zeros and ones (rows)


Aeq=sparse(coordrowseq, coordcolumnseq, fillAeq, numeqconstraints, sU);   %[#equalityconstraints]x[sU]

%%%%%%%%%%%%%%%%%%%%%%  Inequality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Repeat the way we constructed Tcoord_specific but now using the parameter values under consideration
vectors = cellfun(@(x) {x(g,:)}, Ugridreduced); %cell 1xn_U_sets
%For j=1,...,n_U_sets, it picks the g-th row of Ugridreduced{j}
T_temp = cell(1,n_U_sets);
[T_temp{:}] = ndgrid(vectors{:}); 
T_temp = cat(n_U_sets+1, T_temp{:});
T = reshape(T_temp,[],n_U_sets); %matrix of dimension as Tcoord reporting the elements of all the possible n_U_sets-tuples from the U-sets
D=[T(indices_pairs(:,1),:) T(indices_pairs(:,2),:) ...
   Tcoord(indices_pairs(:,1),:) Tcoord(indices_pairs(:,2),:)]; %n_U_sets-tuples n_U_sets-tuples coordinates coordinates
    
%Save in D1 the rows of D_temp where the first triplet is <= then the second
D1_all=D(sum(D(:,1:3)<=D(:,4:6),2)==3, :);
D1=D1_all(:, 7:end); %only coordinates
sd1=size(D1,1);
%Remember <= means: (<,<,<) or (<,=,=) or (<,<,=) or (=,=,=) [any order]
%Identify zero boxes
D1_values= D1_all(:, 1:6);
C1=(round(D1_values(:,1),8)>=round(D1_values(:,5)+D1_values(:,6),8) | round(D1_values(:,4),8)<=round(D1_values(:,2)+D1_values(:,3),8));


%Save in D2 the rows of D_temp where the first triplet is >= (excluding = = =) then the second
D2_all=D((sum(D(:,1:3)>D(:,4:6),2)==3) |... %>,>,>
       (sum(D(:,1:3)>D(:,4:6),2)==2 & sum(D(:,1:3)==D(:,4:6),2)==1) | ... %>,>,=
       (sum(D(:,1:3)>D(:,4:6),2)==1 & sum(D(:,1:3)==D(:,4:6),2)==2), :); %>,=,= 
D2=D2_all(:, 7:end); %only coordinates   
sd2=size(D2,1);
%Identify zero boxes
D2_values=[D2_all(:, 4:6) D2_all(:, 1:3)] ; 
C2=(round(D2_values(:,1),8)>=round(D2_values(:,5)+D2_values(:,6),8) | round(D2_values(:,4),8)<=round(D2_values(:,2)+D2_values(:,3),8));


coordrowsineq=[kron((1:1:sd1), ones(1,8))  kron((sd1+1:1:sd1+sd2), ones(1,8))]; %1x[(sd1+sd2)*8]
%8 comes from the number of all possible vertices generated by each pair of triplets
coordcolumns1=[changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,2)),idU3_t(D1(:,3))) ... 
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,2)),idU3_t(D1(:,3))) ...
                   changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,5)),idU3_t(D1(:,3))) ...
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,5)),idU3_t(D1(:,3))) ...
                   changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,2)),idU3_t(D1(:,6))) ...
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,2)),idU3_t(D1(:,6))) ...
                   changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,5)),idU3_t(D1(:,6))) ...
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,5)),idU3_t(D1(:,6)))]; %sd1x8
coordcolumns2=[changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,2)),idU3_t(D2(:,3))) ... 
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,2)),idU3_t(D2(:,3))) ...
                   changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,5)),idU3_t(D2(:,3))) ...
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,5)),idU3_t(D2(:,3))) ...
                   changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,2)),idU3_t(D2(:,6))) ...
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,2)),idU3_t(D2(:,6))) ...
                   changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,5)),idU3_t(D2(:,6))) ...
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,5)),idU3_t(D2(:,6)))]; %sd2x8
          
coordcolumnsineq=[reshape(coordcolumns1.', 1, 8*sd1) reshape(coordcolumns2.', 1, 8*sd2)];   %transform coordcolumns1,2 in a row vector by reshaping row-wise  
%1x[(sd1+sd2)*8]


fillAineq=[repmat([1 -1 -1 1 -1 1 1 -1], 1, sd1) repmat([-1 1 1 -1 1 -1 -1 1], 1, sd2)]; %remember: signs are flipped because inequalities are <= 
%1x[(sd1+sd2)*8]

bineq=zeros(sd1+sd2, 1);


Aineq=sparse(coordrowsineq, coordcolumnsineq, fillAineq, sd1+sd2, sU);   %[#inequalityconstraints]x[sU]

lower_bineq=-Inf(size(Aineq,1),1);
lower_bineq([C1; C2])=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Run LP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


lower=zeros(sU,1)+eps;
upper=ones(sU,1)-eps;
lower(id_01)=0;
upper(id_01)=1;
% MOSEK
param_MOSEK.MSK_IPAR_LOG = 0; 
prob.blx=lower; %lower bound unknowns
prob.ulx=upper; %upper bound unknowns
prob.c=zeros(sU,1); %objective function
prob.a=[Aeq; Aineq];
prob.blc=[beq; lower_bineq];
prob.buc=[beq; bineq];
[~,result]     = mosekopt('minimize echo(0)',prob, param_MOSEK);
if strcmp(result.sol.itr.prosta, 'PRIMAL_AND_DUAL_FEASIBLE') %if the primal is feasible
   exists(g)=1; 
end

end

%Enlarge results to the original grid of parameter vectors
exists=exists(equalorder);
IdSet = [u01grid(logical(exists),:)  u11grid(logical(exists),:)  u21grid(logical(exists),:) ...
         u02grid(logical(exists),:)  u12grid(logical(exists),:)  u22grid(logical(exists),:)];

